clear; clc; close all;

% 参数设置
a_list = linspace(0.8, 1.1, 100); % a 参数范围
q = 1; % 整数阶
b = 1; P = -2; Q = 2;
N_total = 3000; % 轨迹长度
m_le = 10;      % MLE中临近点搜索最小间隔

% 初始值组（x0, y0）
init_conditions = [
    0.1, 0.5;
   -5.0, 0.5;
    0.1, 0.1;
   -1.0, 0.5
];

colors = ['r', 'b', 'g', 'y'];
legends = {
    '(0.1, 0.5)', '(-5, 0.5)', '(0.1, 0.1)', '(-1, 0.5)'
};

MLE_mat = zeros(length(a_list), size(init_conditions,1));

% 使用parfor并行计算
parfor idx = 1:length(a_list)
    a = a_list(idx);
    mle_vals = zeros(1, size(init_conditions,1)); % 临时存储当前a下的MLE
    for ic = 1:size(init_conditions,1)
        x0 = init_conditions(ic,1);
        y0 = init_conditions(ic,2);
        [x, y] = SCLMM(q, a, b, P, Q, x0, y0, N_total);
        traj = [x; y];
        mle_vals(ic) = MLE_Wolf(traj, m_le);
    end
    MLE_mat(idx, :) = mle_vals; % 并行循环外赋值
end

% 绘图
figure; hold on;
for ic = 1:size(init_conditions,1)
    plot(a_list, MLE_mat(:, ic), 'Color', colors(ic), 'LineWidth', 1.5);
end

xlabel('a', 'FontSize', 14);
ylabel('最大李雅普诺夫指数 (MLE)', 'FontSize', 14);
title('Fig.9(a) 2D-SCLMM 最大李雅普诺夫指数', 'FontSize', 16);
legend(legends, 'Location', 'best');
grid on;
box on;

%% Wolf方法计算MLE示例函数
function MLE = MLE_Wolf(X, m)
    N = size(X, 2);
    MLE_sum = 0;
    count = 0;

    for i = 2:N-m
        dists = sqrt(sum((X(:,i) - X).^2, 1));
        dists(i) = inf;
        idx_min = max(1, i-m);
        idx_max = min(N, i+m);
        dists(idx_min:idx_max) = inf;

        [dist_min, k] = min(dists);

        if dist_min < 1e-12 || isinf(dist_min)
            continue;
        end

        if (i+m <= N) && (k+m <= N)
            dist1 = norm(X(:,i) - X(:,k));
            dist2 = norm(X(:,i+m) - X(:,k+m));
            if dist1 < 1e-12 || dist2 < 1e-12
                continue;
            end
            Ld = log(dist2 / dist1);
            MLE_sum = MLE_sum + Ld;
            count = count + 1;
        end
    end

    if count == 0
        MLE = NaN;
    else
        MLE = MLE_sum / (count * m);
    end
end